Reading data

---

Anteriormente, em MANIPULATING FITS FILES, entendemos como ler os arquivos .fits, providos pelo CoRoT. Todas as curvas de luz com exoplanetas confirmados pelo CoRoT, execeto aquelas que possuiam menos que 584 pontos, forma lidas, reamostradas e salvas em arquivos .csv. Após essa etapa, todos os dados foram disponibilizados online nesse website. Dessa forma, qualquer um pode ter acesso aos dados e importa-los de forma rápida e fácil em qualquer script python.

Além disso, com tantos dados salvos e com a necessidade de consolidar alguns algoritmos em prol de facilitar sua posterior analíse, como ferramentas de plotagem, implementação das filtragens, calcular o phase-fold das curvas de luz, entre outras funciolidades, foi criada uma biblioteca chamada IMT-LightCurve .

Esse Python Package implementa os principais métodos para análise das curvas de luz de uma forma facilitada, amigável e bem documentada, como visto na página de quickstart na sua documentação oficial. A partir desse ponto, utilizaremos suas classes e métodos para auxiliar no projeto.

# !pip install /content/imt_lightcurve-1.2-py3-none-any.whl --force-reinstall
!pip install /content/imt_lightcurve-1.2-py3-none-any.whl

# Importing packages

from imt_lightcurve.models.lightcurve import LightCurve

import pandas as pd
import numpy as np

Utilizando a IMT-LightCurve para carregar os dados disponíveis no GitHub, a leitura das curvas de luz fica padronizada: basta escolher a curva e escrever o nome do arquivo na variável global LIGHTCURVE:str e rodar o bloco abaixo, que já faz o plot da curva escolhida.

# Chosen lightcurve
LIGHTCURVE = 'RESAMPLED_0101086161_20070516T060226'

# Importing lightcurve data from github
data = pd.read_csv('https://raw.githubusercontent.com/Guilherme-SSB/IC-CoRoT_Kepler/main/resampled_files/' + LIGHTCURVE + '.csv')
time = data.DATE.to_numpy()
flux = data.WHITEFLUX.to_numpy()

# Create the LightCurve object
curve = LightCurve(time=time, flux=flux)
print(curve)
curve.plot()

LightCurve Object

O próximo passo do projeto é explorar as diferentes técnicas de filtragens passa-baixa: Filtro de Butterworth, Bessel, Gaussian, Ideal e Median.

Utilizando a descrição do processo de filtragem no domínio da frequência e dos filtros não-lineares, implementaremos cada filtro em funções Python.

Apenas para recordação, a equação que define a filtragem no domínio da frequência é:

\[\begin{split}g(x) = \mathcal{F}^{-1}[H(u)F(u)] \\\end{split}\]

where \(\mathcal{F}^{-1}\) is the Inverse discrete Fourier transform, \(F(u)\) is the Fourier transform of the given function \(f(x)\) (input), \(H(u)\) is the filter transfer function, and \(g(x)\) is the filtered signal (output). Both \(F\), \(H\) and \(G\) are arrays of size \(M\), the same as the input signal. The product \(H(u)F(u)\) is formed using array multiplication; that is, \(G(i, k)=H(i, k)F(i, k)\).

Dessa forma, esse capítulo se baseia em tomar as funções de transferência de cada filtro, discretizá-la e aplicar à modelagem teórica.

Ideal Lowpass Filter - Traduzir

---

Vamos começar com o filtro passa-baixa mais simples de todos, o Filtro Ideal.

Utilizando como base o livro Digital Image Processing, by Gonzalez and Woods, mais precisamente o Capítulo 4.8.1, a definição matemática do filtro é: sendo \(D_0\) uma constante positiva, definida pelo usuário, e \(D(u)\) a distância entre um ponto \(u\) até o centro do espectro de frequência do sinal de entrada, o filtro passa-baixa ideal (ILPF) permite passar todas as componentes de frequência em um círculo de raio \(D_0\) a partir da origem do espectro e remove todas as frequências fora desse círculo. Dessa forma, ele é descrito como:

\[\begin{split} H(u) = \begin{cases} 1, &\text{if } D(u) \le D_0 \\ 0, &\text{if } D(u) \ge D_0 \end{cases} \\\end{split}\]

Como \(D(u)\) é a distância entre um ponto \(u\) até o centro do retângulo de frequência, é definido por

\[\begin{split}D(u) = (u-P/2) \\\end{split}\]

sendo \(P\) o tamanho do vetor original preenchido (padded).

O ponto de transição entre \(H(u) = 1\) e \(H(u) = 0\) é chamado de frequência de corte.

Nota. Na implementação no IMT-LightCurve das funções aqui definidas, será inserida a parcela :

y_filtered += (array.mean() - y_filtered.mean())

visando eliminar um possível problema de deslocamento vertical do resultado da filtragem (que foi observado e acontece em algumas combinações específicas de parâmetros de filtragem, como o filtro de Bessel de ordem 1 e frequência de corte 0.1 Nyquist). Como apenas nos importa os números relativos do fluxo, essa linha de código “alinha” a curva original e a filtrada e melhora a visualização do resultado da filtragem.

def ideal_filter(array, cutoff_freq):
  # Extracting info from curve
  n_time = len(array)
  D0 = cutoff_freq*n_time

  # Procedures to apply the ideal lowpass filter
  expanded = expand_edges(flux, numExpansion=numExpansion)
  fourier = fourier_transform(expanded)

  # Creating the low-pass transfer function array
  i = 0
  for i in range(len(fourier)):
      if fourier[i] > D0:
          fourier[i] = 0

  ifft = inverse_fourier_transform(fourier)
  array_filtered = remove_expand_edges(ifft, numExpansion=numExpansion)
  array_filtered += (flux.mean() - array_filtered.mean())

  return array_filtered

Nota. Como o Filtro passa-baixa Ideal não requereu função de transferência, o retorno da função já é o array de entrada filtrado. As demais técnicas, que exigem função de transferência, irão retornar o array do filtro, que ainda deve sofrer todos os procedimentos descritos no Capítulo 2 - Filters.

Assim como nas demais implementações, as funções expand_edges, fourier_transform, inverse_fourier_transform e remove_expand_edges foram explicadas e implementadas no junto a descrição teórica dos processos de filtragem

Plotting some results

Vamos variar a frequência de corte do filtro e verificar o resultado.

Aqui vale uma explicação também: a classe LightCurve do nosso Python Package criado, possui a implementação de todas as técnicas aqui descritas. Na documentação há uma explicação detalhada de como foi feita essa implementação e como utilizar todos os recursos da biblioteca .Dessa forma, a aplicação dos procedimentos se tornou tão simples quanto visto abaixo:

filtered = curve.ideal_lowpass_filter(cutoff_freq=0.1)
filtered.view_filtering_results()

Gaussian Lowpass Filter

---

O próximo filtro que estudaremos é o Filtro Gaussiano passa-baixas. Por mais que não tenha sido necessário o calculo da função de transferência do Filtro passa-baixa Ideal, o usual é utilizar a função de transferência discreta dos filtros. A partir de agora, todas as técnicas no domínio da frequência (Gaussian, Butterworth and Bessel Low-pass Filters) requerão a discretização de sua função de transferência e, para isso, continuaremos a utilizar como base as definições do livro do Gonzalez and Woods.

Dessa forma, the transfer function of a Gaussian 1-D lowpass filter (GLPFs) is defined by

\[H(u) = e^{-D^2(u)/2D^2_0}\]

where \(D(u)\) was defined on Equation (1) and \(D_0\) is a positive constant, which can be interpreted as the cutoff frequency.

def gaussian_array(array, cutoff_freq):
  # Extracting info from curve
  n_time = len(flux)
  D0 = cutoff_freq*n_time
  xc = n_time

  # Procedures to filtering on frequency domain
  expanded = expand_edges(flux, numExpansion=numExpansion)
  padded = padding(expanded)
  centralized = centralize_fourier(padded)
  fourier = fourier_transform(centralized)

  # Creating the low-pass transfer function array
  len_filter = len(fourier)
  filter_array = np.zeros(len_filter)

  i = 0
  for i in range(len_filter):
      filter_array[i] = exp((-(i-(xc-1.0))**2)/(2*((D0)**2)))

  return filter_array

Note. The cutoff frequency must be given in Nyquist.

Plotting some results

filtered = curve.gaussian_lowpass_filter(cutoff_freq=0.1)
filtered.view_filter_results()

filtered = curve.gaussian_lowpass_filter(cutoff_freq=0.4)
filtered.view_filter_results()

Butterworth Lowpass Filter

---

The transfer function of a Butterworth 1-D lowpass filter (BLPF) of order \(n\), and with cutoff frequency at a distance \(D_0\) from the origin, is defined as

\[H(u) = \frac{1}{ 1 + [D(u) / D_0]^{2n} }\]

where \(D(u)\) and \(D_0\) was defined on Eq. (1).

By the definition, the Butterworth filter have two free parameters: the cutoff frequency and the filtering order. Then, we can modify both, as we can see on the code cell below, intending to have the best results possibles.

Note. The cutoff frequency must be given in Nyquist.

def butterworth_array(array, fourier_transform, cutoff_freq, order):
  # Extrating information of the signal
  n_time = len(array)
  D0 = cutoff_freq * n_time
  xc = n_time

  # Creating the filter array
  len_filter = len(fourier_transform)
  filter = np.zeros(len_filter)

  for i in range(len_filter):
    filter[i] = 1.0 / (1.0+(abs(i-(xc-1.0))/D0)**(2.0*order))

  return filter

Plotting some results

filtered = curve.butterworth_lowpass_filter(2, 0.1)
filtered.view_filter_results()

filtered = curve.butterworth_lowpass_filter(4, 0.3)
filtered.view_filter_results()

Bessel Lowpass Filter

---

The transfer function, \(H(s)\), of a Bessel lowpass filter is defined by

\[\begin{split}H(s) = \frac{\theta_n (0) }{\theta_n (s/ \omega_0)} \\\end{split}\]

where

\[\theta_n (s) = \sum_{k=0}^{n} a_k s^k\]

and

\[\begin{split}a_k = \frac{(2n-k)!}{2^{n-k}k!(n-k)!} \\\end{split}\]

Showing how it works…

Parameters

order = 2
cutoff_freq = 0.6

Control lib

from control import *

### Computing ak

from math import factorial

coef = []
i = 0
while i <= order:
  ak = (factorial(2*order - i)) / ( 2**(order - i)*factorial(i)*factorial(order - i) )
  # print(ak)
  coef.append(ak)
  i += 1

print(coef)

[3.0, 3.0, 1.0]

### Computing θn(s)

s = TransferFunction.s
theta_array = []
k = 0
for k in range(order+1):
  theta_n = coef[k] * (s**k)
  theta_array.append(theta_n)
  # numerical_numerator = coef[0]
  # print(theta_n)

print(theta_array[0])
print(theta_array[1])
print(theta_array[2])

3
-
1


3 s
---
 1


s^2
---
 1

### Computing H(s)

coef_numerator = theta_array[0]

list_denominator = theta_array[:]

denominator = 0
for item in list_denominator:
  denominator += item

print(denominator)

s^2 + 3 s + 3
-------------
      1

### Filling in transfer function

G = coef_numerator / denominator

print(G)
print(type(G))

      3
-------------
s^2 + 3 s + 3

<class 'control.xferfcn.TransferFunction'>

Applying

def bessel(array, fourier_transform, cutoff_freq, order):

  # Extracting features from signal
  n_time = len(array)
  D0 = cutoff_freq * n_time
  xc = n_time

  # Creating the bessel transfer function array
  len_filter = len(fourier_transform)
  filter = np.zeros(len_filter)
  i=0

  for i in range(len_filter):
    filter[i] = np.real(evalfr(G, ( np.abs(i-(xc-1.0))/D0 )))

  return filter

Plotting some results

filtered = curve.bessel_lowpass_filter(2, 0.1, numExpansion=100)
filtered.view_filter_results()

filtered = curve.bessel_lowpass_filter(4, 0.3, numExpansion=100)
filtered.view_filter_results()

Median Filter

---

O filtro mediano, por sua vez, é aplicado de uma forma consideravelmente diferente, na qual cada valor dos dados filtrados corresponde a uma mediana de um grupo de valores adjacentes nos dados originais. Esse filtro já é utilizado com certa frequência em análises de curvas de luz.

from scipy.signal import medfilt

def median_filter(array, window_size):
    return medfilt(array, window_size)

Reference: Scipy Documentation

Plotting some results

filtered = curve.median_filter(3)
filtered.view_filter_results()

filtered = curve.median_filter(9)
filtered.view_filter_results()

Saving data

Agora que todos os processos foram detalhados e implementados, vamos aplicar cada técnica de filtragem e diversas combinações de seus respectivos parâmetros em todo o dataset de curvas de luz reamostradas com exoplanetas confirmados.

• Para os casos que têm como parâmetro livre a frequência de corte, será avaliado as frequências:

cutoff_freq = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]

• Para os casos que têm a ordem de filtragem como parâmetro livre, será avaliado os valores:

orders = [1, 2, 3, 4, 5, 6]

• Para o único caso que depende do número de vizinhos, Filtro Mediano, será condirado:

numNeighbors = [3, 5, 7, 9, 11]

• Para os filtros que exigem mais de um parâmetro, será avaliado todas as combinações dos intervalos definidos.

Para salvar os dados, vamos utilizar o método export_filters_to_csv, que pertence à classe LightCurve, basta passar como parâmetro qual técnica de filtragem deseja salvar, seus respectivos parâmetros, onde deseja salvar e onde está o dataset das curvas reamostradas.

curve = LightCurve(time=None, flux=None)

# curve.export_filters_to_csv(WHERE_TO_SAVE, DATASET_PATH, 'ideal', cutoff_freq_range=(0.1, 0.9, 0.1))
# curve.export_filters_to_csv(WHERE_TO_SAVE, DATASET_PATH, 'gaussian', cutoff_freq_range=(0.1, 0.9, 0.1))
# curve.export_filters_to_csv(WHERE_TO_SAVE, DATASET_PATH, 'butterworth', cutoff_freq_range=(0.1, 0.9, 0.1), order_range=(1, 6, 1))
# curve.export_filters_to_csv(WHERE_TO_SAVE, DATASET_PATH, 'bessel', cutoff_freq_range=(0.1, 0.9, 0.1), order_range=(1, 6, 1))
# curve.export_filters_to_csv(WHERE_TO_SAVE, DATASET_PATH, 'median', numNei_range=(3, 11, 2))

API

---

Here, you can apply all the filtering processes, with whatever parameters you want, on whatever lightcurve

Select the Lightcurve:

RESAMPLED_0100725706_20070516T060226

RESAMPLED_0101086161_20070516T060226

RESAMPLED_0101206560_20070516T060226

RESAMPLED_0101368192_20070516T060050

RESAMPLED_0102671819_20071023T223035

RESAMPLED_0102671819_20120112T183055

RESAMPLED_0102708694_20071023T223035

RESAMPLED_0102708694_20120112T183055

RESAMPLED_0102725122_20071023T223035

RESAMPLED_0102725122_20120112T183055

RESAMPLED_0102764809_20071023T223035

RESAMPLED_0102890318_20070206T133547

RESAMPLED_0102912369_20070203T130553

RESAMPLED_0105118236_20100708T204534

RESAMPLED_0105209106_20080415T231048

RESAMPLED_0105228856_20100408T223049

RESAMPLED_0105793995_20080415T231048

RESAMPLED_0105819653_20080415T231048

RESAMPLED_0105833549_20080415T231048

RESAMPLED_0105891283_20080415T231048

RESAMPLED_0106017681_20080415T231048

RESAMPLED_0110839339_20081116T190224

RESAMPLED_0110864907_20081116T190224

RESAMPLED_0221686194_20081011T143035

RESAMPLED_0300001097_20081116T190224

RESAMPLED_0310247220_20090403T220030

RESAMPLED_0311519570_20090403T220030

RESAMPLED_0315198039_20100305T001525

RESAMPLED_0315211361_20100305T001525

RESAMPLED_0315239728_20100305T001525

RESAMPLED_0630831435_20110708T151253

RESAMPLED_0652180928_20110708T151253

RESAMPLED_0652180991_20110708T151253

LIGHTCURVE = 'RESAMPLED_0110839339_20081116T190224'

Filtering processes availables:

• Ideal (Cutoff frequency)

• Gaussian (Cutoff frequency)

• Butterworth (Order, Cutoff frequency)

• Bessel (Order, Cutoff frequency, numExpansion=100)

• Median (Number of neighbors)

from utils import *
import pandas as pd
import numpy as np

Loading BokehJS ...

data = pd.read_csv('https://raw.githubusercontent.com/Guilherme-SSB/IC-CoRoT_Kepler/main/resampled_files/' + LIGHTCURVE + '.csv')
time = data.DATE.to_numpy()
flux = data.WHITEFLUX.to_numpy()

curve = lightcurve.LightCurve(time, flux)
curve.plot()

help(curve.how_to_filter)

Help on method how_to_filter in module utils.lightcurve:

how_to_filter(order: int, cutoff_freq: float, numNei: int, numExpansion: int) method of utils.lightcurve.LightCurve instance
    This function describes how to filtering using this library

    Parameters
    ----------
    order : int
        Used in Butterworth and Bessel filtering. Matches the filter order.

    cutoff_freq : float
        Used in Ideal, Gaussian, Butterworth and Bessel filtering. Matches the cutoff frequency.

    numNei : int
        Used in Median. Matches the number of neighbors to consider

    numExpansion : int
        Used in all processes. Corresponds to how much you want to expanded the curve's edges
        (to avoid some problems caused by the Fast Fourier Transform algorithm).
        Preliminary tests show that all processes works fine with numExpansion=70,
        except for Bessel filtering which required numExpansion=100


    Methods
    -------
    curve.ideal_lowpass_filter(cutoff_freq)

    curve.gaussian_lowpass_filter(cutoff_freq)

    curve.butterworth_lowpass_filter(order, cutoff_freq)

    curve.bessel_lowpass_filter(order, cutoff_freq)

    curve.median_filter(numNei)


    Examples
    --------
    >>> from utils import *
    >>> curve = lightcurve.Lightcurve(time=[1, 2, 3, 4], flux=[10, 15, 12, 20])
    >>> print(curve)
    <utils.lightcurve.LightCurve object>
    >>> curve.plot()

    >>> filtered = curve.bessel_lowpass_filter(order=3, cutoff_freq=0.4, numExpansion=100)
    >>> print(filtered)
    <utils.lightcurve.FilteredLightCurve object>

    >>> filtered.view_filter_results()

filtered = curve.butterworth_lowpass_filter(order=2, cutoff_freq=0.3)
filtered.view_filter_results()